.libPaths()
# delete data in global environment:
rm(list = ls())
# load standard packages:
library("stringr", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("reshape", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("ez", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("ggplot2", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("plyr", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("nortest", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("car", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("ltm", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("lsr", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("afex", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("nlme", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("ggthemes", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("robustHD", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("cowplot", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("tidyr", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("lattice", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")                                
library("plyr", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("WRS2", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("multinomRob", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("vcd", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("pscl", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("Hmisc", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("doBy", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("multcomp", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("arm", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
library("sjPlot", lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")

# load some necessary functions, for within-subject, repeated measure designs
load("/Users/XXX/effect.size.Rdata")
load("/Users/XXX/mean.sd.Rdata")
load("/Users/XXX/min.mean.sd.max.Rdata")
load("/Users/XXX/summarySE.Rdata")
load("/Users/XXX/summarySEwithin.Rdata")
load("/Users/XXX/normDataWithin.RData")
load("/Users/XXX/rcontrast.RData")

#read csv file:
getwd()
setwd("XXX")
data_for_r = read.csv("data_time_exp_future_banana.csv")

#delete NAs
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}
data_for_r<-completeFun(data_for_r, "Response")

# response time: subtract 666ms (= stimulus length)
data_for_r$rt_offset<-data_for_r$rt-666 #  total amount: 6840 - 6814 - 6750
# only Responses that were given before stimulus onset:
data_for_r<-subset(data_for_r, rt_offset>0) # 26 deleted
data_for_r<-subset(data_for_r, rt_offset<6000) # 64 deleted
# 90 out of 6840 outliers: 1.3 %

#select only 'test blocks' (without learning blocks) per condition
data_for_r$test<-grepl(pattern = "(test)", data_for_r$blockcode) 
data_for_r<-subset(data_for_r, test=="TRUE")  

# rename conditions:
data_for_r<-mutate(data_for_r,blockcode=as.character(blockcode))
data_for_r<-mutate(data_for_r,blockcode=sapply(strsplit(data_for_r$blockcode, split='_', fixed=TRUE),function(x) (x[2])))
data_for_r$blockcode=gsub(data_for_r$blockcode, pattern = "vorwaerts", replacement = "forward")
data_for_r$blockcode=gsub(data_for_r$blockcode, pattern = "rueckwaerts", replacement = "backward")
data_for_r$blockcode=gsub(data_for_r$blockcode, pattern = "stehen", replacement = "standing")
data_for_r$condition <-data_for_r$blockcode

# only trials were an answer (correct or incorrect) was given
data_correct_and_incorrect<-subset(data_for_r, correct=="correct" | correct=="incorrect")  
# only correct answers:
data_correct<-subset(data_for_r, correct=="correct") 

#select relevant columns:
data_correct<-data_correct[c("participant", "Response", "condition", "correct", "order")]
data_correct_and_incorrect<-data_correct_and_incorrect[c("participant", "Response", "condition", "correct", "order")]

# count number of correct and incorrect items:
data_correct_count<-count(data_correct, c("correct", "participant", "condition", "Response", "order"))
data_correct_and_incorrect_count<-count(data_correct_and_incorrect, c("correct", "participant", "condition", "Response", "order"))

# density plots per Response & condition
ggplot(data_correct_count, aes(x = freq)) + geom_density() + facet_wrap(Response ~ condition)

# descriptive
datac<-summaryBy(freq ~ condition*Response, data = data_correct_count, 
                 FUN = function(x) { c(m = mean(x), s = sd(x), var = var(x)) } )
datac # -> no obvious differences

# function, atop, numbers
external_graphs = function(ext = TRUE){
  if( is.rstudio() ){
    if(isTRUE(ext)){
      o = tolower(Sys.info()["sysname"])
      a = switch(o,
                 "darwin"  = "quartz",
                 "linux"   = "x11",
                 "windows" = "windows")
      options("device" = a)
    } else{
      options("device"="RStudioGD")
    }
    
    # Kill open graphic devices
    graphics.off()
  }
}


ggplot(datac, aes(x=condition, y=freq.m, fill=Response)) +
  geom_bar(position=position_dodge(.9), colour="black", stat="identity") +
  geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=freq.m-freq.s, ymax=freq.m+freq.s)) +
  theme_classic() +
  theme(legend.position="top") +
  xlab("\nCondition") +
  ylab(expression(atop("Amount of correct items (mean)", atop(italic("Possible range of scores: 0 - 10"), "")))) +
  labs("d") +
  theme(text = element_text(size=50))+
  coord_cartesian(ylim=c(0,10)) +
  scale_fill_manual(values=c("#CCCCCC","#FFFFFF")) +
  theme_bw() +
  theme(text = element_text(size=30))+
  geom_hline(yintercept=38) 


################## model: backward vs standing, forward vs standing #######################

# We create a model where we compare each condition against the control condition
data_correct_and_incorrect$condition <- factor(data_correct_and_incorrect$condition, levels = c("standing", "backward", "forward"))
data_correct_and_incorrect$condition<-as.factor(data_correct_and_incorrect$condition)
contrasts(data_correct_and_incorrect$condition)

# We conduct a generalized linear mixed model with a binomial distribution 
# see also http://ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/mixed%20model%20guide.html
# baseline model:
summary(m_baseline <- glmer(correct ~ 1 + (1 | participant) + (1+condition | participant) + (1+Response | participant), data = data_correct_and_incorrect,  family = binomial(link="logit"))) 

# adding response
summary(m_Response <- glmer(correct ~ Response + (1 | participant) + (1+condition | participant) + (1+Response | participant), data = data_correct_and_incorrect,  family = binomial(link="logit"))) 
# compare baseline and response model
anova(m_baseline, m_Response)
# -> Adding Response does not improve the model.

# adding condition:
summary(m_condition <- glmer(correct ~ condition + (1 | participant) + (1+condition | participant) + (1+Response | participant), data = data_correct_and_incorrect,  family = binomial(link="logit"))) 
# compare baseline and condition model
anova(m_baseline, m_condition)
# descriptive
(data_condition<-summaryBy(freq ~ condition, data = data_correct_count, 
                 FUN = function(x) { c(m = mean(x), s = sd(x), var = var(x)) } ))
# -> Main effect condition: Participants in the bakcward condition had slightly less correct answers (mean = 7.03) compared to the control condition (mean = 7.38). 

# model with response + condition
summary(m_response_and_condition <- glmer(correct ~ Response+condition + (1 | participant) + (1+condition | participant) + (1+Response | participant), data = data_correct_and_incorrect,  family = binomial(link="logit"))) 
# model with interaction response * condition
summary(m_response_and_condition_interaction <- glmer(correct ~ Response*condition + (1 | participant) + (1+condition | participant) + (1+Response | participant), data = data_correct_and_incorrect,  family = binomial(link="logit"))) 
anova(m_response_and_condition, m_response_and_condition_interaction)
# -> The interaction term does not improve the model

# check fixed effects correlation matrix
sjp.glmer(m_response_and_condition_interaction, type = "fe.cor")

# check qqplot of random effects
sjp.glmer(m_response_and_condition_interaction, type = "re.qq")

# As expected from the descriptive statistics: 
# m_response_and_condition_interaction: No main effects of condition, Response, and no interaction effects.

